/* 
 * File:   NumerovMatrixSparse.cpp
 * Author: murmeli
 * 
 * Created on March 12, 2013, 7:51 AM
 */

#include "NumerovMatrixSparse.h"

NumerovMatrixSparse::NumerovMatrixSparse() {
    setName("Numerov matrix sparse GPU");
    setMatrixMethod(true);
    setSparse(true);
}

NumerovMatrixSparse::NumerovMatrixSparse(const NumerovMatrixSparse& orig) {
}

NumerovMatrixSparse::~NumerovMatrixSparse() {
}

void NumerovMatrixSparse::methodStep(double& step, vector<double>& xx, vector<double>& yy, vector<double>& dy, cDY* DY, cPot* V) {
    int elements = (int)(V->getUniverse()/step);
    double VStart = V->getPotLeft();
    double sq = step*step;
    vector<double> tempx;
    //Eigen::LDLT<Eigen::MatrixXd> solver;
    Eigen::SimplicialLDLT< Eigen::SparseMatrix<double> > solver;
    Eigen::SparseMatrix<double> mA(elements,elements),mB(elements,elements),mV(elements,elements), temp(elements,elements);
    
   // mA.setZero();
   // mB.setZero();
   // mV.setZero();
   // temp.setZero();
 
    tempx.reserve(elements);
    
    tripletsA.reserve(elements);
    tripletsB.reserve(elements);
    tripletsV.reserve(elements);

    tripletsA.push_back(T(0,0,-2.0/sq));
    tripletsA.push_back(T(0,1,1.0/sq));
    
    tripletsB.push_back(T(0,0,10.0/12));
    tripletsB.push_back(T(0,1,1.0/12));
    
    tempx.push_back(V->getPotLeft());    
    
    tripletsV.push_back(T(0,0,V->getPotentialValueAt(VStart)));
    
    //cout << elements << endl;
    
    for(int i = 1; i < elements-1; i++ ) {
        tripletsA.push_back(T(i,i-1,1.0/sq));
        tripletsA.push_back(T(i,i,-2.0/sq));
        tripletsA.push_back(T(i,i+1,1.0/sq));
        
        tripletsB.push_back(T(i,i-1,1.0/12.0));
        tripletsB.push_back(T(i,i,10.0/12.0));
        tripletsB.push_back(T(i,i+1,1.0/12.0));
        
        tripletsV.push_back(T(i,i,V->getPotentialValueAt(VStart+i*step)));
//        cout << i << endl;
    tempx.push_back(tempx.back() + step);  
}

    int i = elements - 1;
    tripletsA.push_back(T(i,i-1,1.0/sq));
    tripletsA.push_back(T(i,i,-2.0/sq));
        
    tripletsB.push_back(T(i,i-1,1.0/12.0));
    tripletsB.push_back(T(i,i,10.0/12.0));
    tripletsV.push_back(T(i,i,V->getPotentialValueAt(VStart+i*step)));

    tempx.push_back(tempx.back() + step);
    
    mA.setFromTriplets(tripletsA.begin(), tripletsA.end());
    mB.setFromTriplets(tripletsB.begin(), tripletsB.end());
    mV.setFromTriplets(tripletsV.begin(), tripletsV.end());
    
    tripletsA.clear();
    tripletsB.clear();
    tripletsV.clear();

    solver.compute(mB);
    temp = -0.5 * solver.solve(mA);
    temp+=mV;
    setSparseH(temp);
    //setH(temp);
    setX(tempx);
}